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ABSTRACT 

We determine the velocity vector of M31 with respect to the Milky Way 
and use this to constrain the mass of the Local Group, based on Hubble Space 
Telescope proper-motion measurements of three fields presented in Paper I. We 
construct A-body models for M31 to correct the measurements for the contri- 
butions from stellar motions internal to M31. This yields an unbiased estimate 
for the M31 center-of-mass motion. We also estimate the center-of-mass mo- 
tion independently, using the kinematics of satellite galaxies of M31 and the 
Local Group, following previous work but with an expanded satellite sample. 
All estimates are mutually consistent, and imply a weighted average M31 helio- 
centric transverse velocity of (vw,vn) = (—125.2 ± 30.8, —73.8 ± 28.4) kms -1 . 
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We correct for the reflex motion of the Sun using the most recent insights into 
the solar motion within the Milky Way, which imply a larger azimuthal veloc- 
ity than previously believed. This implies a radial velocity of M31 with respect 
to the Milky Way of Kad,M3i — —109.3 ± 4.4 kms" 1 , and a tangential velocity 
Kan,M3i = 17.0 kms -1 , with la confidence region Vtan,M3i < 34.3 kms -1 . Hence, 
the velocity vector of M31 is statistically consistent with a radial (head-on col- 
lision) orbit towards the Milky Way. We revise prior estimates for the Local 
Group timing mass, including corrections for cosmic bias and scatter, and obtain 
M LG = M M w,vir + M M 3i, v ir = (4.93 ± 1.63) x 10 12 M . Summing known esti- 
mates for the individual masses of M31 and the Milky Way obtained from other 
dynamical methods yields smaller uncertainties. Bayesian combination of the 
different estimates demonstrates that the timing argument has too much (cos- 
mic) scatter to help much in reducing uncertainties on the Local Group mass, 
but its inclusion does tend to increase other estimates by ~ 10%. We derive a 
final estimate for the Local Group mass from literature and new considerations 
of M LG = (3.17 ± 0.57) x 10 12 M . The velocity and mass results imply at 95% 
confidence that M33 is bound to M31, consistent with expectation from observed 
tidal deformations. 

Subject headings: galaxies: kinematics and dynamics — Local Group — M31. 

1. Introduction 

The Milky Way (MW) is a member of a small group of galaxies called the Local Group 
(LG). The LG is dominated by its two largest galaxies, the MW and the Andromeda galaxy 
(M31). The mass and dynamics of this group have been the topic of many previous studies 
(e.g., van den Bergh 2000; van der Marel & Guhathakurta 2008, hereafter vdMG08; Li & 
White 2008; Cox & Loeb 2008; and references therein). Analysis of these topics is important 
for interpretation of structures inside the LG, such as dark halos, satellite galaxies, and tidal 
streams. It is also important for understanding the LG in a proper cosmo logical context, 
since it provides the nearest example of both Large Scale Structure and hierarchical galaxy 
formation. While much progress has been made in understanding the LG mass and dynamics, 
this has not been based on actual knowledge of the three-dimensional velocity vector of M31. 
This is because until now, the proper motion (PM) of M31 has been too small to measure 
with available techniques. 

In Paper I (Sohn, Anderson & van der Marel 2012) we reported the very first absolute 
PMs of M31 stars in three different fields observed with the Hubble Space Telescope (HST): 
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a field along the minor axis sampling primarily the M31 spheroid (the "spheroid field"), a 
field along the major axis sampling primarily the M31 outer disk (the "disk field"), and a 
field along M31's Giant Southern Stream (GSS) sampling primarily the stars that constitute 
this stream (the "stream field"). For each field we measured the average PM of the M31 stars 
with respect to the stationary reference frame of background galaxies. The results are listed 
in Table [TJ PMs in mas/yr were converted to velocities (vw,vn) in km/s in the directions 
of West and North using the known distance D of M31. Throughout this paper we adopt 
D = 770 ± 40 kpc (see references in vdMG08). The velocity uncertainties are dominated by 
the PM uncertainties, with distance uncertainties making only a minimal contribution. 

In the present paper we use the observed PMs to determine both the direction and 
size of the M31 velocity vector with respect to the MW, and we use this knowledge with 
the Local Group timing argument (Kahn & Woltjer 1959; vdMG08; and Li & White 2008) 
to estimate the LG mass. We then compare the velocity and mass results to independent 
estimates of the same quantities. For example, vdMG08 estimated the transverse motion of 
M31 based on the kinematics of satellite galaxies of M31 and the Local Group. Furthermore, 
the mass of the Local Group has been estimated independently by adding up the individual 
masses of the MW and M31, as estimated from various dynamical tracers (e.g. Klypin et 
al. 2002; Watkins et al. 2010). By statistically combining all the results we are able to build 
an improved and comprehensive understanding of the dynamics and mass of the LG. 

The outline of this paper is as follows. In Section 2 we use iV-body models of M31 and its 
prominent tidal substructures to calculate predictions for the internal kinematics of M31 stars 
in the three fields observed with HST. We use the results to correct the transverse velocities 
measured with HST, to estimate the transverse velocity of the M31 center-of-mass (COM). 
In Section 3 we revisit the methods of vdMG08 to estimate the M31 transverse motion from 
the kinematics of satellites, but with an expanded satellite sample. We combine the results 
with the HST measurements to obtain a final estimate for the M31 transverse motion. In 
Section 4 we derive the corresponding space motion in the Galactocentric rest frame, taking 
into account the latest insights about the solar motion in the MW. The results are consistent 
with a radial orbit for M31 towards the MW. In Section 5 we use the M31 motion to estimate 
the LG mass using the timing argument. We find that the estimate is quite uncertain due 
to cosmic scatter, and we show how a more accurate estimate can be obtained by combining 
it with estimates of the individual MW and M31 masses. In Section 6 we consider the 
galaxy M33, the third most massive galaxy of the Local Group (van den Bergh 2000), and 
we derive its relative velocity with respect to M31. We also derive an estimate for the mass 
of M33, and show that M33 is most likely bound to M31, as is usually assumed. We use 
this knowledge to further refine our estimate for the Local Group mass. In Section 7 we 
discuss and summarize the main results of the paper. An Appendix presents a discussion of 
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various parameterizations used in the literature (and the paper text) to quantify the dark 
halo density profiles and masses of galaxies. Where necessary to compare the properties 
of Local Group galaxies with predictions from cosmological simulations, we use a Hubble 
constant H = 70kms~ 1 Mpc -1 and a matter density Vt m = 0.27 (Jarosik et al. 2011). 

This is the second paper in a series of three. Paper III (van der Marel et al. 2012, in 
prep.) will present a study of the future orbital evolution and merging of M31, M33, and 
the MW, using the velocities and masses derived in the present paper as starting conditions. 

2. Correction for Internal Kinematics 

The PMs measured with HST in M31 fields contain contributions from both the M31 
COM motion, and from the internal kinematics of M31. In each field, different fractions 
of the stars are contributed by different structural components. Specifically, the galaxy has 
different equilibrium components, including both a disk and spheroids (bulge/halo). We will 
refer to these jointly as the "base galaxy". The galaxy also contains material that is in the 
process of being accreted. This includes in particular the material responsible for the creation 
of the GSS (which in fact is spread out over a large fraction of the projected area of the 
galaxy, and not just the actual position of the Stream). To estimate the M31 COM motion, 
we need to know for each field observed with HST both the fraction of the stars in each 
component, and the transverse motion kinematics of those stars. The fractional contributions 
can in principle be estimated purely observationally from line-of-sight (LOS) velocity studies. 
However, estimates of the transverse motion kinematics requires a full dynamical model, since 
these motions are not directly constrained observationally. We therefore resort to N-body 
models for M31 like those previously constructed by some of us (e.g., Geehan et al. 2006; 
Fardal et al. 2006, 2007, 2008) to understand various observed features of M31. 

2.1. N-body models of M31 Structure 

The M31 model we use here is constructed from two separate but related parts. The 
base galaxy is an N-body realization of a model of the mass and light in M31 itself. The 
GSS component is a snapshot from a dynamical N-body simulation of the formation of the 
GSS, performed using the same mass model of M31. Taken together, these two components 
reproduce reasonably well the features in M31 that are expected to contribute to our HST 
fields. 

The base galaxy, which is a slightly altered version of the model from Geehan et 



-5 - 



al. (2006), contains bulge, disk, and halo components. The bulge and disk are assumed 
to be mostly baryonic and therefore trace the light. To the dark halo present in Geehan 
et al. (2006) we add a stellar halo, which is necessary to reproduce the extended power-law 
component that has been discovered in the halo regions (Guhathakurta et al. 2005; Irwin 
et al. 2005; Kalirai et al. 2006b; Chapman et al. 2006; Ibata et al. 2007). We assume this 
stellar halo follows the mass distribution of the dark halo, although it contains only a tiny 
fraction of that mass. When added together, these components satisfy the surface-brightness 
profiles of M31's bulge, disk, and halo regions reasonably well. Most importantly for this 
study, they also satisfy a series of kinematic constraints, including the disk rotation curve, 
the bulge velocity dispersion, and constraints on the halo mass from statistical tracers such 
as globular clusters, planetary nebulae, satellite galaxies, and red giant stars. We created 
the particle realization of this model using the ZENO library (Barnes 2011). 

The GSS component is created by simulating the disruption of a satellite galaxy, in a 
fixed potential corresponding to the mass model just discussed. The model is an updated 
version of that found in Fardal et al. (2007), to which we refer the reader for a physical 
discussion. This model uses a spherical progenitor, although it is possible that the progenitor 
may in fact have been a disk galaxy (Fardal et al. 2008). After starting at large radius with 
carefully chosen initial conditions, the satellite disrupts at its first pericentric passage. The 
model is evolved using the PKDGRAV tree N-body code (Stadel 2001) for nearly 1 Gyr, 
until it forms orbital wraps closely resembling features in M31 including the GSS itself, and 
the NE and W shelves. We refer to all the particles generated by this component as the GSS 
component, regardless of where they are on the sky. 

All the parameters of the base model, N-body simulations, and data-model comparisons 
are presented in Fardal et al. (2012, in prep.). However, many properties and details are 
similar to preceding papers (Geehan et al. 2006; Fardal et al. 2006, 2007, 2008). Besides the 
morphological evidence, the GSS model satisfies a set of observational constraints, includ- 
ing the detailed kinematic pattern in the W Shelf, the precise sky position of the GSS, the 
distance to various fields along the GSS (McConnachie et al. 2003), and their peak LOS ve- 
locities (Ibata et al. 2004; Guhathakurta et al. 2006; Kalirai et al. 2006a; Gilbert et al. 2009). 
We do not use the observed color-magnitude diagrams (CMDs) of the HST PM fields to con- 
strain the models, since those are not easily decomposed into distinct structural components. 
In fact, the CMDs of the HST spheroid field and the HST stream field are strikingly similar, 
given that they are believed to be dominated by different structural components (Brown et 
al. 2006). 

Figure [1] shows a smoothed projected view of the N-body model. The GSS is visible 
South-East of the galaxy center, and the North-East and Western shelf are emphasized with 
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dashed outlines. This image can be compared to star-count maps of giant stars in M31, 
which show the same features (e.g., Ibata et al. 2005; Gilbert et al. 2009; or Paper I, which 
also shows the location of our PM fields). Figure |5] shows a smoothed view of the N-body 
model in LOS velocity vs. projected distance space, for particles South-East of the galaxy 
center. There is good agreement between the outline of the GSS in this representation (dark 
band in the figure), and the observed peak LOS velocity of the GSS as a function of radius 
(blue points), including that measured in the HST stream PM field of Paper I (circle). 

2.2. Proper-Motion Corrections for HST Fields 

The predictions of the N-body model for M31 are summarized in Table [1] for each of 
the three HST fields. The quantity /base is the fraction of the stars that belongs to the base 
galaxy, and fcss is the fraction of the stars that belongs to the GSS. The average velocities 
in the LOS, W, and N directions are listed for both the base and the GSS components, and 
also for their properly weighted average ("all"). The quantities were extracted over fields 
that are somewhat larger than the HST fields (up to 10 arcmin from the field center), to 
decrease the N-body shot noise. All velocities are expressed in a reference frame in which 
M31 is at rest. 

The average internal transverse velocity kinematics of the M31 stars in the HST fields 
(v\y, Vn){oII) are generally small, always below 125 km/s in absolute value. There are several 
reasons for this. In the HST spheroid field we are sampling primarily the spheroidal compo- 
nents of M31. At large radii these have limited mean rotation (Dorman et al. 2012). In the 
HST disk field we are primarily sampling the M31 disk, which has a sizeable circular velocity 
(~ 250 km/s). However, M31 has a large inclination. So along the major axis, most of the 
rotation is seen along the LOS, and not in the transverse direction. In the HST stream field 
we are primarily sampling the GSS, which has a significant mean three-dimensional velocity 
(254 km/s). However, the inclination of the stream is such that most of the velocity is seen 
along the LOS. Moreover, some 20% of the stars in the stream field do not belong to the 
GSS, but mostly to the spheroid component. 

For all three fields, the contribution to the observed transverse motion from the internal 
kinematics of M31 is similar to or smaller than the random errors in the HST measurements. 
Even significant fractional changes in the model predictions therefore do not strongly affect 
our final results. It therefore did not seem worthwhile in the context of the present study 
to further refine the model. Nonetheless, it is worth pointing out that the model is far 
from perfect, and that there are some salient features of M31 that it fails to reproduce. 
For example, the model with a spherical progenitor overestimates the contribution of GSS 
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particles on their first orbital wrap to fields along the minor axis (especially those more 
distant than our spheroid field). Also, LOS velocity studies of the GSS have revealed a 
secondary cold component (in addition to the GSS and base galaxy; e.g., Kalirai et al. 2006a; 
Gilbert et al. 2007, 2009) which is not reproduced by our model. This could be, e.g., from 
a severely warped disk component, or from wrapped-around material of a GSS loop not 
included in our model. And finally, some authors have proposed models for the structure 
of M31's outer and accreted components that differ from those in our models (e.g., Ibata et 
al. 2005; Gilbert et al. 2007). 

To estimate the transverse M31 COM motion from the data for each HST field, we 
first subtract the contribution from internal M31 kinematics (v(all)) from the measurement 
(v(HST)). We then correct for the effect of viewing perspective as described in van der 
Marel et al. (2002) and vdM08. This corrects for the fact that at the position of each field, 
the M31 COM systemic LOS and transverse velocity, respectively, are not exactly aligned 
along the local LOS and transverse directions. The corrections are small (below 10 km/s 
in absolute value), because all fields are located within 2 degrees of the M31 center. The 
final estimates are listed in Tabled] as (vw, vn)(COM), and are summarized also at the top 
of Table [3J In propagating the uncertainties, we assigned an uncertainty of 20 km/s per 
coordinate to the model of the M31 internal kinematics. This number need not be known 
accurately, since the final uncertainties are always dominated by measurement errors in the 
PMs. 

The results for the three different fields are mutually consistent with each other at the 
la level (see also Paper 10). This justifies the use of a straightforward weighted average to 
combine the results, which gives (v w ,v N )(COM) = (-162.8 ± 47.0, -117.2 ± 45.0) km/s. 
For comparison, the direct weighted average of the HST observations, with no corrections for 
internal M31 kinematics, is (v w ,v N )(HST) = (-154.1 ±44.9, -112.9±42.9) km/s. Clearly, 
the corrections for internal kinematics make only a small difference for the final transverse 
velocity estimate. The fact that the differences are below 10 km/s is due to our combination 
of results for well-chosen fields, since the per-field corrections are much larger than this. 

3. Transverse Velocity from Satellite Kinematics 

In vdMG08 we presented several methods for estimating the space motion of M31 from 
the kinematics of satellites, which assume that the satellites of M31 and the LG on average 



1 Paper I defines a x 2 quantity, ^| , that assesses the extent to which measurements for different fields 
agree. Table [3] implies x| = 3.5, for TVdf = 4 degrees of freedom. 
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follow their motion through space. The M31 transverse velocity derived in that paper has 
random error bars of 34 to 41 km/s. This is somewhat smaller than what we have obtained 
here from the HST PM measurements, although the systematic error bars on the vdMG08 
values may be larger (because the underlying methodology makes more assumptions). Either 
way, these results remain of considerable interest as an independent constraint on the M31 
space motion. We therefore update here the results from vdMG08 using additional satellite 
data that has become available more recently. 

3.1. Constraints from Line-of-Sight Velocities of M31 Satellites 

The first method of vdM08 is based on the fact that any transverse motion of M31 
induces an apparent solid body rotation in the line-of-sight velocity field of its satellites, 
superposed on their otherwise primarily random motions. The amplitude and major axis 
of the rotation field are determined by the size and direction of M31's transverse motion. 
In vdM08 we constrained the M31 transverse motion by fitting the velocities of 17 M31 
satellites with known line-of-sight velocities. 

For the present study we added the satellites listed in Table [2j These are objects that 
previously either did not have LOS velocity measurements available, or which had not yet 
been discovered. This includes six dSph galaxies: And XI, XIII, XV, XVI, XXI, and XXII. 
Three other recently discovered dSph galaxies, And XVII, XIX, XX, have not yet had their 
LOS velocities measured. As in vdMG08, we leave out And XII and XIV, because their 
large negative LOS velocities with respect to M31 indicate that they may be falling into 
M31 for the first time (Chapman et al. 2007; Majewski et al. 2007). We also leave out the 
more recently discovered And XVIII, which may be too distant from M31 to be directly 
associated with it (McConnachie et al. 2008). We do include And XXII, even though it 
may be a satellite of M33 rather than M31 (Martin et al. 2009; Tollerud et al. 2012). We 
note that And IV is not included in our combined sample because it is a background galaxy 
(Ferguson et al. 2000), while And VIII is not a galaxy at all (Merrett et al. 2006). For all 
dSphs, including those from Paper I and not listed in Table [2j we used the newly measured 
LOS velocities from Tollerud et al. (2012), where available. Otherwise, the values listed in 
Paper I or the sources listed in Table |2] were used. Our new sample in Table [2] also includes 
the 8 globular clusters of M31 that lie at projected distances > 40 kpc and have known LOS 
velocities. 

We repeated the vdMG08 analysis, using the combined sample of the 17 satellites in 
Table 1 of vdMG08 and the 14 satellites in Table [2j The implied space motion of M31 is 
listed in TableOin the row labeled "M31 satellites". The result for (vw, vn) differs from that 
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derived in vdMG08 by (—40.1, 13.4) km/s. This is considerably smaller than the error bars 
in the result of (144.1,85.4) km/s. The addition of the 14 new satellites has not decreased 
the error bars on the result. This is in part because most satellites are observed relatively 
close to M31, so that any solid-body rotation signal would be small. As before, the new 
result is roughly consistent at the la level with zero transverse motion. So no solid-body 
rotation component is confidently detected, which in turn implies that M31 cannot have a 
very large transverse motion. The fits imply a one-dimensional velocity dispersion for the 
satellite sample of <j sat = 84.8 ± 10.6 km/s. This is 8.5 kms -1 larger than the value derived 
in vdMG08, which again is within the uncertainties. 

3.2. Constraints from Proper Motions of M31 Satellites 

The second method of vdM08 is based on the M31 satellites M33 and IC 10. These 
galaxies have accurately known PMs from water-maser observations (Brunthaler et al. 2005, 
2007). The three-dimensional velocity vectors of these galaxies give an estimate of the M31 
space motion to within an accuracy of <T sat per coordinate. Transplanting the M33 and IC 
10 velocity vectors to the position of M31, followed by projection onto the local LOS, W, 
and N directions, yields the results listed in Table |3] in the rows labeled "M33 PM" and 
"IC 10 PM". These are identical to what was derived in vdMG08, but with slightly larger 
uncertainties (due to the increased estimate of cr sat ). 

3.3. Constraints from Line-of-Sight Velocities of Outer Local Group Galaxies 

The third method of vdMG08 is based on the line-of-sight velocities of Local Group 
satellites that are not individually bound to the MW or M31. In vdMG08 the method 
was applied to 5 satellites (see their Table 2). The Cetus dSph (RA= 6.54597°, DEC= 
—11.04432°) was excluded because of lack of knowledge of its LOS velocity at the time. For 
the present study we have rerun the analysis including Cetus, using t>Los = — 87 ± 2 km/s 
from Lewis et al. (2007). Its distance D = 755 ± 23 kpc (McConnachie et al. 2005) places 
Cetus at -D bary ps 600 kpc from the Local group barycenter. With addition of the Cetus dSph 
to the vdMG08 analysis, the implied space motion of M31 is listed in Table [3] in the row 
labeled "Outer LG Galaxies". The result for (vw,vn) differs from that derived in vdMG08 
by (—14.9,-13.5) km/s. This is considerably smaller than the error bars in the result of 
(58.0,52.5) km/s. 
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3.4. Comparison and Combination of Constraints 

In general, modeling of satellite galaxy dynamics can be complicated for a variety of 
reasons, especially when the goal is to estimate galaxy masses: the satellite system may not 
be virialized, with continueing orbital evolution (Mateo et al. 2008) or infall (e.g., Chapman 
et al. 2007; Majewski et al. 2007); the distribution of satellite orbits may not be isotropic 
(e.g., Watkins et al. 2010); satellites on large-period orbits are not expected to be randomly 
distributed in orbital phase (e.g., Zaritsky & White 1994); satellites may have correlated 
kinematics (e.g., van den Bergh 1998); and the three-dimensional distribution of satellites 
may not be spherical (Koch & Grebel 2006) or symmetric (McConnachie & Irwin 2006). 
However, many of these potential issues do not affect the analysis that we have presented here 
and in vdMG08 to estimate the M31 transverse velocity. Sections 13.11 and 13.21 only assume 
that the M31 satellites are drawn from a distribution that has the same mean velocity as 
M31, and which has no mean rotation. Section [3731 only assumes that the LG satellites are 
drawn from a distribution that has the same mean velocity as the LG barycenter. Virialized 
equilibrium, isotropy, random phases, or symmetry are not required. Nonetheless, there 
is always the possibility that residual systematics may have affected the results. To get a 
handle on this, we have compared in detail the results for the M31 transverse velocity from 
the different techniques. 

The {vw,vn) for M31 inferred from the different methods in Sections I3.1H3.3I are in 
mutual agreement to within the uncertainties. The same was true also in the original analysis 
of vdMG08. Since the methods and the underlying data are quite different for the various 
estimates, this in itself is a direct indication that any residual systematics cannot be large. 
Since the results from the different methods are in agreement, it is reasonable to take their 
weighted average, as listed in Table [3J The result for (vw,vn) differs from that derived in 
vdMG08 by (—19.0,-7.1) km/s. This is considerably smaller than the error bars in the 
result of (40.7,36.6) km/s, so the new analysis presented here has not significantly altered 
the results previously derived by vdMG08. 

An even stronger check on any residual systematics is provided by Figure |3J It compares 
the weighted average of the HST PM measurements (with corrections for internal kinematics) 
from Section [2] (as listed in Table |3]) with the weighted average from the updated vdMG08 
analysis. The difference between these results is (Avw, Av^) = (—65.8 ±62.2, —72.1 ±58.0). 
This means that the results are consistent within the uncertainties: the probability of a 
residual this large occurring by chance in a two-dimensional Gaussian distribution is 26%. 
Since the methods employed are totally different, and have quite different scopes for possible 
systematic errors, this is very successful agreement. This suggests not only that there are 
no large residual systematics in the results from the satellite kinematics, but also that there 



- 11 - 



are no large residual systematics in the M31 PM analysis. This is a very important cross- 
check, since the displacements on which our PM measurements are based are below 0.01 
detector pixels (see Paper I for a detailed discussion of the systematic error control in the 
PM analysis). 

Since the HST PM analysis and the satellite kinematics analysis yield statistically con- 
sistent results for the M31 transverse velocity, it is reasonable to take the weighted average 
of the results from the two methodologies. This yields 

(%,%) = (-125.2 ±30.8, -73.8 ± 28.4) kms" 1 , (1) 

as listed in the bottom row of Table E] and shown in black in Figure [3J This is the final result 
that we use for the remainder of our analysis. 

4. Space Motion 

4.1. Galactocentric Rest Frame and Solar Motion 

As in vdMG08, we adopt a Cartesian coordinate system (X,Y,Z), with the origin 
at the Galactic Center, the Z-axis pointing towards the Galactic North Pole, the X-axis 
pointing in the direction from the Sun to the Galactic Center, and the F-axis pointing in the 
direction of the Sun's Galactic rotation. We choose the origin of the frame to be at rest (the 
Galactocentric rest frame), and we wish to calculate the velocity of galaxies in this frame. 
This requires knowledge of the solar velocity in the Milky Way, since the solar reflex motion 
contributes to any observed velocities (such as the heliocentric velocities listed in Table |3]). 

In vdMG08 we adopted the standard IAU values (Kerr & Lynden-Bell 1986) for the 
distance of the Sun from the Galactic Center Rq = 8.5 kpc, and the circular velocity of the 
Local Standard of Rest (LSR), Vq = 220 km/s. Neither of these quantities has historically 
been known particularly accurately though, and their exact values continue to be debated. 
Recently, a number of new methodologies have become available. These provide new insights 
into Ro and Vo, and we therefore use the results from these studies here. 

Some of the best constraints on Rq now come from studies of the orbits of stars around 
the Sgr A* supermassive black hole. Gillessen et al. (2009) obtained Ro = 8.33 ± 0.35 kpc 
(consistent also with Ghez et al. 2008). Most of the available constraints on the velocity Vo 
are actually constraints on the ratio Vq/Rq. The best constraint on this ratio now comes 
from the observed PM of Sgr A*, since the black hole is believed to be at rest in the 
galaxy to within ~ 1 km/s. Reid & Brunthaler (2004) obtained that (Vo + Vp ec )/-Ro = 
30.2 ± 0.2km s~ 1 kpc _1 . Here V pec is the peculiar velocity of the Sun in the rotation direction. 
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In vdMG08 we adopted the solar peculiar velocity from Dehnen & Binney (1998). However, 
there is now increasing evidence that V pec from that study (and other studies) is too small 
by ~ 7 km/s. We adopt here the more recent estimates from Schonrich, Binney, & Dehnen 
(2010): (U pec , V pec , W pec ) = (11.1,12.24,7.25), with uncertainties of (1.23, 2.05, 0.62) kms" 1 
(being the quadrature sum of the random and systematic errors). Combination of these 
results implies that Vq = 239.3 ± 10.3 km/s, significantly larger than the canonical IAU 
value of 220 km/s. The uncertainty is dominated entirely by the uncertainty in R , and the 
errors in Vq and Rq are highly correlated. 

Observations of masers in high-mass star-formation regions in the MW have been used to 
argue independently for a value of Vq in excess of the canonical 220 km s" 1 (Reid et al. 2009). 
However, McMillan & Binney (2010) showed that these data by themselves do not strongly 
constrain the Galactic parameters. On the other hand, McMillan (2011) showed that when 
combined with the other constraints described above through detailed models, the maser data 
do help to constrain Rq more tightly, and therefore Vq. He obtained: Rq = 8.29 ± 0.16 kpc 
and Vq = 239 ± 5 km/s. These are the values we adopt here. 



4.2. M31 Space Motion 

Based on the adopted M31 distance and solar parameters, the position of M31 in the 
Galactocentric rest frame is 

f M31 = (-378.9, 612.7, -283.1) kpc. (2) 

The velocity of the Sun projects to (v sys , vw, vn)q = (191.9, 142.5, 78.5)kms _1 at the position 
of M31. Since one observes the reflex of this, these values must be added to the observed 
M31 velocities to obtain its velocity in the Galactocentric rest frame. The velocity vector 
corresponding to the observed COM LOS velocity t>Los = — 301 ± 1 kms -1 (vdMG08) and 
the final weighted average (vw,vn) given in Tabled transformed to the Galactocentric rest 
frame, is then 

v M31 = (66.1 ±26.7, -76.3 ± 19.0, 45.1 ±26.5) kms" 1 . (3) 

The errors (which are correlated between the different components) were obtained by prop- 
agation of the errors in the individual position and velocity quantities (including those for 
the Sun) using a Monte-Carlo scheme. 

If the transverse velocity of M31 in the Galactocentric rest frame, V tan , equals zero, then 
M31 moves straight towards the Milky Way on a purely radial (head-on collision) orbit. This 
orbit has (fvK,^A r )rad = (—141.5 ± 3.0, —78.8 ± 1.7) kms -1 (this is approximately the reflex 
of the velocity of the Sun quoted above, because the lines from the Sun to M31 and from the 
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Galactic Center to M31 are almost parallel). The listed uncertainty is due to propagation 
of the uncertainties in the solar velocity vector. The radial orbit is indicated as a starred 
symbol in Figure [3J The velocity VM31 calculated in the previous paragraph corresponds to 
a total velocity 

\V mi \ = 110.6 ±7.8 kms -1 . (4) 

The radial velocity component is 

Kad,M3i = -109.2 ±4.4 km s- 1 , (5) 

and the tangential velocity component is 

Vtan,M3i = 17.0 kms -1 (la confidence region : Vt an ,M3i < 34.3 kms" 1 ). (6) 

The uncertainties were calculated as in vdMG08, using a flat Bayesian prior probability for 
V iaa . These results imply that the velocity of M31 is statistically consistent with a radial 
orbit at the la levelj^] 

It has been known for a long time that the transverse velocity of M31 is probably 
Vtan,M3i ^ 200km s _1 . The large scale structure outside the LG does not provide enough tidal 
torque to have induced a much larger transverse motion, subsequent to the radial expansion 
started by the Big Bang (Gott & Thuan 1978; Raychaudhury & Lynden-Bell 1989). Li & 
White (2008) used the ACDM cosmological Millenium simulation to identify some thousand 
galaxy pairs that resemble the MW-M31 pair in terms of morphology, isolation, circular 
velocities, and radial approach velocity. The median tangential velocity of the pairs was 
Kan — 86km s _1 , with 24% of the pairs having V tlin < 50km s _1 (see their figure 6). Therefore, 
the observed M31 tangential velocity is somewhat below average compared to cosmological 
expectation, but it is not unusually low. The exact reason why the tangential velocity of the 
MW-M31 pair has ended up below-average is not clear, but it may be related to the details 
of the growth history and the local environment of the LG. 

Peebles et al. (2001) showed that at values Vt a n,M3i ^ 200km s _1 , many velocities can be 
consistent with the observed positions and velocities of galaxies in the nearby Universe. Our 
new observational result that Vt an ,M3i < 34.3 kms -1 at la confidence therefore significantly 
reduces the parameter space of possible orbits. Peebles et al. (2011) recently proposed 
a model for the history and dynamics of the LG in which Vt an ,M3i — 100.1 kms -1 and 
(vw,vn) = (—240.5, —63.1) kms -1 . This is inconsistent with our final velocity estimates at 
> 3a confidence. 



2 If the IAU value Vq = 220 km/s is used for the LSR circular velocity and the Dehnen & Binney 
(1998) values are used for the solar peculiar velocity (as in vdMG08), then a radial orbit has (%, t>Ar) ra d = 
(—126.6, —71.4) kms -1 . With these assumptions, the inferred velocity of M31 is still statistically consistent 
with a radial orbit at the la level. 
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5. Local Group Mass 

The velocity vector of M31 with respect to the Milky Way constrains the mass of the 
Local Group through the so-called "timing argument" (Kahn & Woltjer 1959; Lynden-Bell 
1981, 1999; Einasto & Lynden-Bell 1982; Sandage 1986; Raychaudhury & Lynden-Bell 1989; 
Kroeker & Carlberg 1991; Kochanek 1996). Recent applications of this method were pre- 
sented in vdMG08 and Li & White (2008). In Section [5.11 we provide a revised estimate 
of the Local Group timing mass using the new insights into the M31 velocity vector from 
Section 14.21 In Section 15.21 we combine the result in statistical fashion with results from 
other independent methods for constraining the Local Group mass. 

Different studies often quote different mass quantities. However, for proper use and 
comparison it is important to transform all measurements to a common definition. For each 
galaxy, the total mass is dominated by a very extended dark halo. Common characterizations 
of dark matter halos are summarized in Appendix |A] The density profile is often modeled as 
an NFW profile (Navarro et al. 1997; eq. |A2j ) or a Hernquist (1990; eq. |A8j ) profile. The 
former has infinite mass, while the latter has finite mass Mh- Common characterizations of 
halo masses also include the virial mass M vir enclosed within radius r vir (eq. |Alj ). the mass 
M 2 oo enclosed within the radius r 2 oo (eq. |A5j ). or the mass M(r) enclosed within some given 
physical radius r in kpc (eqs. |A3|IA9j ). 

In our discussion of galaxy masses, we transform all results into M v j r estimates. The 
transformation requires knowledge of the density profile, which for this purpose we assume 
to be of the NFW form with known concentration c vir (eq. |A3j ). For the MW and M31 we 
take c v j r = 10 ± 2, based on a combination of specific models (Klypin et al. 2002; Besla et 
al. 2007), and cosmological simulation results (Neto et al. 2007; Klypin et al. 2011). This 
implies M 200 /M viv = 0.839 ± 0.014 (eq. \A7\). 

5.1. Timing Argument 

Under the assumption of Keplerian motion, the relative orbit of M31 and the MW 
is determined by four parameters: the total mass M tot ; the semi- major axis length a (or 
alternatively, the orbital period T); the orbital eccentricity e; and the current position within 
the orbit, determined by the eccentric anomaly 77. In turn, four observables are available to 
constrain the orbit: the current M31 distance D; the radial and tangential M31 velocities in 
the Galactocentric rest frame, K a d,M3i and V tlint M3i] and the time t since the last pericenter 
passage, which should be equal to the age of the Universe t = 13. 75 ±0.11 Gyr (Jarosik et 
al. 2011; when the matter of both galaxies originated together in the Big Bang). There are 
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as many observables as unknowns, so the orbital parameters and Mlg can be determined 
uniquely. The implied value of M tot is called the "timing mass" . Most commonly, the relevant 
equations are solved under the assumption of a radial orbit (e = 1 and VJ; ani M3i = 0), but 
this assumption is not necessary when a measurement of Kan,M3i is actually available (e.g., 
vdMG08). 

For the M31 space motion derived above, the "timing mass" M tot timing = (4.27 ±0.53) x 
1O 12 M . When a radial orbit is assumed (which is consistent with the data) then M tot timing = 
(4.23 ± 0.45) x 10 12 M (this is somewhat smaller, because any transverse motion increases 
the timing mass). The listed uncertainties are the RMS scatter of Monte-Carlo simulations 
as in vdMG08, which take into account all the observational uncertainties. The black curve 
in the top panel of Figure H] shows the complete probability histogram for the radial orbit 
case. 

These timing mass results are about 1.0-1.3 x 10 12 M Q lower than what was obtained 
by, e.g., vdMG08 and Li & White (2008). This can be viewed as an improvement, since 
previous estimates of M tot appeared anomalously high compared to independent estimates 
of the masses of the individual M31 and MW galaxies (as summarized in vdMG08). Some 
of the decrease in timing mass is due to the fact that Vt an ,M3i found here is slightly smaller 
than in vdMG08. However, most of the decrease in timing mass is not due to the new HST 
measurements, but due to the new values for the solar motion used here. The solar velocity 
in the F-direction, vy = V + V pec , is 251.2 km/s in our calculations here. By contrast, it was 
26 km/s lower in the calculations of vdMG08. The Y component of the solar motion projects 
predominantly along the LOS direction towards M31, and not the W and N directions. The 
component of the solar motion in the LOS direction is therefore 20.2 km/s higher than what 
it was in vdMG08. As a consequence, in the Galactocentric rest frame, M31 approaches 
the MW with a radial velocity that is 20.2 km/s slower than what is was in vdMG08. This 
slower approach implies a lower timing mass. 

The timing argument equations are based on a simple Keplerian formalism. To assess 
how accurate this argument is in a cosmological context one must make comparisons to N- 
body simulations (Kroeker & Carlberg 1991). Li & White (2008) did this for the currently 
favored ACDM cosmology, using results from the Millennium simulation. They identified 
simulated galaxy pairs like the MW-M31 system, with known masses, and quantified the 
accuracy of the radial orbit timing argument. They found that the timing argument has 
very little bias, when viewed as an estimate of the sum M tot 20 o of the galaxy's M 2 oo values, 
but significant scatter. They quantified this "cosmic scatter" , by averaging over pairs with 
all possible transverse velocities. However, not surprisingly, their Figure 6 shows that the 
scatter increase with V tan . Since we now know that M31 actually has a low value of V ta , n , it 
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is more appropriate to quantify the cosmic scatter by restricting the statistics to pairs in the 
simulation with low V tim . We measured by hand from their Figure 6 all pairs with V taun < 50 
km/s, and extracted the ratio M totj 2oo/^tot,timin g - We folded the probability distribution of 
these ratios into our Monte-Carlo scheme for estimating the total mass from the radial orbit 
timing argument. This yields the estimate Mt t,2oo — (4.14 ± 1.36) x 10 12 M , which has a 
three times larger uncertainty than what is implied by observational errors alone. This can 
be converted into an estimate for the summed virial masses using to formulae of Appendix A, 
which yields 

M tot) vir = (4.93 ± 1.63) x 10 12 M (timing argument). (7) 

This is our final estimate from the timing argument, which takes into account observational 
uncertainties, cosmic bias, and cosmic scatter. The red curve in the top panel of Figure 0] 
shows the complete probability histogram for M tot vir . 



5.2. Combination with Other Milky Way and M31 Mass Constraints 

The best alternative method for estimating the mass of the Local Group is to add up 
estimates of the masses of the individual M31 and MW galaxies^] Estimates for the masses 
of these galaxies were already summarized in vdMG08, so here we highlight primarily some 
more recent results. 

Watkins et al. (2010) studied the kinematics of M31 satellites and found that the mass 
within 300 kpc is determined fairly robustly, M M3 i(300kpc) = (1.40 ±0.43) x 10 12 M Q . The 
quoted uncertainty is the quadrature sum of the random error of 0.40 x 10 12 M , and a sys- 
tematic uncertainty of 0.15 x 1O 12 M due to the assumed velocity anisotropy of the satellites. 
There may be other systematic uncertainties in the analysis, but these are more difficult to 
quantify and are neglected here. At this mass and with the relevant halo concentrations, 
M(300kpc)/M vir = 1.018±0.002 (see Appendix[Aj). Hence, M M3 i,vir = (1.38±0.43) x 1O 12 M . 
We use this to set a Gaussian probability distribution for M M 3 liVir in our discussion below 
(blue curve in top panel of Figure H]). This is consistent with the study of Klypin et al. (2002), 
which folded in a wider range of observational constraints, and obtained successful models 
with dark halo masses M M 3i,dark,vir of either 1.43 x 10 12 M or either 1.60 x 10 12 M . This 
corresponds to a total mass of M M3 i,vir = 1-52 x 10 12 M or 1.69 x 10 12 M , respectively, 
after adding in also the combined stellar mass of the M31 disk and bulge. The Watkins et 



3 The other method of estimating M to t from the size of the Local Group turn-around radius yields esti- 
mates that tend to be biased low. A radial infall model is generally assumed, which is almost certainly an 
oversimplification (see vdMG08). 
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al. results are also consistent with the recent results of Tollerud et al. (2012). They applied 
a mass estimator calibrated on cosmo logical simulations to the M31 satellite kinematics and 
obtained M M3 i, v ir = l-2±g;? x 10 12 M Q . 

Watkins et al. (2010) showed that mass estimates for the mass of the MW from satellite 
kinematics are much more uncertain. This is due to the unknown velocity anisotropy, com- 
bined with the fact that we see most satellites almost radially from near the Galactic Center. 
Good mass estimates therefore need to fold in a more diverse set of observational constraints. 
Moreover, uncertainties are reduced significantly by assuming that the radial profile of the 
dark matter is known, and follows a cosmologically motivated parameterization. McMillan 
(2011) used such methods to obtain a dark halo mass MMw,dark,2oo = (1.26±0.24) x 1O 12 M . 
This corresponds to M MW dark vir = (1.50 ±0.29) x 1O 12 M . This is consistent with the study 
of Klypin et al. (2002), who favored M M w,dark,vir = 1-0 x 1O 12 M , but showed that reasonable 
models with MMw,dark,vir = 2.0 x 10 12 M can be constructed as well. Adding the combined 
stellar mass of the MW disk and bulge to obtain the total Mmw.vu- adds ~ 0.06 x 10 12 to 
these values. The rapid motion of the Magellanic Clouds and Leo I have been used to argue 
for masses at the high end of this range of values (e.g., Zaritsky et al. 1989; Shattow & 
Loeb 2008; Li & White 2008; Boylan-Kolchin et al. 2011). However, the underlying assump- 
tions in these arguments cause significant uncertainties. Based on the range of results in 
the literature, we adopt here, fairly arbitrarily, a flat probability distribution for Mmw.vu- 
between 0.75 and 2.25 x 10 12 M (green curve in top panel of Figure H]). This distribution 
has the same mean (1.50 x 1O 12 M ) as inferred by McMillan (2011), and the same dispersion 
(0.43 x 10 12 Mq) as we use for M31, but with a broader, flatter shaped 

We use the listed probability distributions for the individual M31 and MW masses 
as priors; this also sets a prior probability distribution for M tot vir = M M 3i, v ir + Mmw,™ 
(magenta curve in top panel of Figure H]). We then fold in the timing argument results 
to determine posterior probability distributions, as follows. We draw a random mass from 
the probability distribution for Mt t,vir derived from the timing argument (red curve in top 
panel of Figure H]). We then draw a random MM\y,vir from its prior distribution. We then 
calculate the corresponding M M3l!Vir = M tot vir — M MW vir , and its probability p given the 
prior distribution for M M 3i iV ir- This set of values is then accepted or rejected in Monte-Carlo 



4 There do exist models in the literature that yield or use higher mass estimates for M31 and the MW 
than we use here. This includes estimates based on halo occupation distributions (e.g., Guo et al. 2010) or 
the timing argument (e.g., Loeb et al. 2005; Cox & Loeb 2008). It should be kept in mind though that such 
estimates are statistical in nature. Cosmic scatter must therefore be taken into account, and this yields large 
uncertainties (e.g., Li & White 2008; Guo et al. 2010). It is therefore important that any mass estimate for 
an individual galaxy, as opposed to an ensemble of galaxies, also take into account the actually observed 
resolved properties, rotation curves, and satellite kinematics. 
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sense, depending on the probability p. We thus build up posterior probability distributions 
for M M w,vir, Mvi3i,vir, and M tot vir , which are shown as dotted lines with the same colors in 
the bottom panel of Figure HI 

The posterior distribution for MM3i, V ir is still roughly Gaussian, but its average has 
increased from (1.38 ±0.43) x 1O 12 M to (1.51 ±0.42) x 1O 12 M . The posterior distribution 
for M MW vir is not flat like its prior, but skewed towards higher masses. Its average has 
increased from 1.50 to 1.63 x 10 12 M . The likelihood of MW and M31 masses at the low 
end of the prior distributions is significantly reduced in the posterior distributions. This is 
relevant, since some previously reported mass estimates do fall on this low- mass end (e.g., 
Evans et al. 2000; Ibata et al. 2004). At the best-estimate virial masses for the MW and 
M31, the corresponding virial radii are 308 kpc and 300 kpc, respectively (eq. |Alj ). Since 
the distance between the galaxies is D = 770 ± 40 kpc, the virial spheres are not currently 
overlapping. 

The prior distribution for M totjV i r corresponds to (2.88 ± 0.61) x 10 12 M , while its 
posterior distribution corresponds to (3.14 ± 0.58) x 10 12 M . Therefore, inclusion of the 
timing argument increases the estimate of the LG mass by only ~ 9%, due to the large 
cosmic variance. Since this is considerably smaller than the prior uncertainties on Mmw.vu- 
and MM3i,vir, the timing argument does not in fact help much to constrain the total LG 
mass, beyond what we already know from the MW and M31 individually. We have found 
this to be a robust conclusion, independent of the exact probability distributions adopted 
for M31 and the MW, and independent of the exact solar and M31 motion adopted in the 
timing argument. 



6. Mass Constraints from M33 

The galaxy M33 is the most massive companion of M31 (e.g., van den Bergh 2000). 
In the past decade, evidence has been found from both HI (Braun & Thilker 2004) and 
star-count maps (McConnachie et al. 2009) for tidal features indicative of past interactions 
between these galaxies. Models for these features such as those presented by McConnachie et 
al. (2009) require that M33 be bound to M31@ The galaxy M33 is one of the few galaxies in 
the Local Group for which an accurate PM measurement is available from VLBA observations 
of water masers. Hence, combined with our new M31 results, the relative motion of M33 
with respect to M31 is now known with reasonable accuracy (Section l6.ip . The mass of M33 



5 Whilc the MW also has a massive companion, namely the Large Magellanic Cloud, it is unclear whether 
the galaxies in this pair form a bound system (Besla et al. 2007). 
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can also be estimated independently (Section 16.2 j) . With knowledge of the relative velocity 
and mass, the assumption that M33 is bound to M31 can be used to further refine our 
understanding of the M31 mass, and hence the Local Group mass (Section 16. 3p . 

6.1. M33 Space Motion 

To establish the binding energy of the M31-M33 system, we need to know the cur- 
rent position r*M33 and velocity VM33 of M33 in the Galactocentric rest frame. These were 
determined in similar fashion as for M31 (see Section |4T2"|) . but now based on the following 
observables: a distance D M33 = 794±23kpc (McConnachie et al. 2004), line-of-sight velocity 
^ios,M33 = —180 ± 1 kins -1 (vdMG08), and PM from water masers as measured by Brun- 
thaler et al. (2005) and discussed in vdMG08. This yields f M3 3 = (-476.1, 491.1, -412.9)kpc, 
and Vu33 = (43.1 ± 21.3, 101.3 ± 23.5, 138.8 ± 28.1) kms" 1 . The observational errors in the 
Galactocentric velocity of M33 are similar to those for M31 reported in Section 14.21 

The positions of the three galaxies MW, M31 and M33 define a plane in the Galacto- 
centric rest frame. For simplicity, we will refer to this plane as the "trigalaxy plane". It is of 
interest for understanding the orbital evolution of the MW-M31-M33 system, to know how 
the M33 velocity vector is oriented with respect to this plane. To assess this, we introduce a 
new Cartesian coordinate system (X', Y', Z') based on the following definitions: the frame 
has the same origin as the (X,Y,Z) system (i.e., the Galactic Center); the X'-axis points 
from the origin to M31 at t = 0; the V'-axis is perpendicular to the X'-axis, and points from 
M31 to M33 as seen in projection from the Galactic Center; and the Z'-axis is perpendicular 
to the X'- and y'-axes in a righthanded sense. With these definitions, the trigalaxy plane is 
the (X', Y') plane. Hence, let us refer to the (X', Y', Z') system as the "trigalaxy coordinate 
system" . 

Based on the position vectors r M31 and rM33 from Sections 14.21 and 16. 1[ the unit vectors 
of the (X', Y r , Z') system can be expressed in (X, Y, Z) coordinates as 

u x , = (-0.48958,0.79153,-0.36577), 
u Y , = (-0.47945, -0.60013, -0.64029), 

u z , =(-0.72632,-0.13810,0.67331). (8) 

If r is a vector expressed in Galactocentric (X, Y, Z) coordinates, then the corresponding 
vector r* expressed in the trigalaxy coordinate system is 

f'=(X',Y',Z') = (f-u x ,,r-u YI ,f-u z ,), (9) 
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where ■ denotes the vector inner product. We use equations fl8|9|) as the fixed definition of 
the (X', Y', Z') system throughout this paper, even when we vary the positions of the Sun, 
M31, and M33 within their observational uncertainties. 

Observational uncertainties of ~ 20-30km s _1 aside, the velocities of M31 and M33 in tri- 
galaxy coordinates are = (—109.2, —15.5, —7. l)kms _1 and-i^m — (8.3, —170.3, 48. 2)kms 
These vectors make angles with the (X',Y') plane of only —3.7° and 15.8°, respectively. By 
definition, the MW galaxy currently has zero velocity in the Galactocentric rest frame. How- 
ever, the gravitational attraction from M31 and M33 will set it in motion with a velocity 
directed in the (X', Y') plane. Hence, all three galaxies start out in the (X', Y') plane, 
with velocity vectors that are close to this plane. This implies that the orbital evolution of 
the entire MW-M31-M33 system will happen close to the trigalaxy plane, with the "verti- 
cal" Z'-component playing only a secondary role. Detailed calculations of the future orbital 
evolution and merging of the MW-M31-M33 system are the topic of Paper III. 



6.2. M33 Mass 

The mass of M33 is not negligible with respect to that of M31. It is therefore necessary 
to know the mass of M33 to determine whether the M31-M33 system is bound. Corbelli 
(2003) modeled the rotation curve and mass content of M33. The rotation curve rises to 
~ 130 km s _1 out to the last data point at 15 kpc. Since the data do not reveal a turnover 
in the rotation curve, both the halo concentration and virial mass are poorly constrained 
(Fig. 6b of Corbelli 2003). Moreover, the rotation field is complex with significant twisting 
(Corbelli & Schneider 1997). This complicates interpretation in terms of circular motion. 
To estimate the M33 virial mass it is therefore necessary to use more indirect arguments. 
For this, we compare M33 to M31. 

Corbelli (2003) used her rotation curve fits to estimate the mass-to-light ratio of the 
M33 disk. From this, she inferred a stellar mass 2.8-5.6 xlO 9 M at 3a confidence^ Higher 
values correspond to a maximum-disk fit, while lower values correspond to a sub-maximal 
disk. Guo et al. (2010) instead used the observed B — V color of M33 with stellar population 
model predictions to estimate the mass-to-light ratio. With an assigned uncertainty of 0.1 
dex for this method, one obtains (2.84 ± 0.73) x 10 9 M . We combine these methods into a 
single rough estimate M M3 3,* = (3.2 ± 0.4) x 10 9 M . 



6 Corbelli's mass scale for Ho = 65kms -1 Mpc -1 was transformed to the Hubble constant used here. The 
small mass contribution from the nuclear component of M33 is well within the quoted uncertainties. M33 
has no bulge. 
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For M31, Klypin et al. (2002) used rot at ion- curve fits to estimate both the disk and the 
bulge mass. The two models they present cover the ranges M M31disk = 7.0-9.0 x 10 10 M 
and M M31bulgc = 1.9-2.4 x 10 10 M . Upon subtraction of the gas mass of ~ 0.6 x 10 10 (van 
den Bergh 2000), this yields M M3 i,* = (8.3-10.8) x 10 10 M . The Guo et al. (2010) method 
based on the galaxy B — V color yields instead Mm3i,* = (7.0 ± 1.8) x 1O 1O M . We combine 
these methods into a single rough estimate Mm3i,* = (7.9 ± 0.9) x 10 10 M . 

These estimates imply that Mm33,*/Mm3i,* = 0.041 ± 0.007. This can be compared 
to the baryonic mass ratio implied by the Tully Fisher relation, (VM33/VM31) 4 (McGaugh 
2005). With V M 33 « 130 km s" 1 (Corbelli & Salucci 2000) and V M3 i ~ 250 km s" 1 (Corbelli 
et al. 2010) this yields 0.073. This is consistent with the estimate of Mm33,*/Mm3i,*, if one 
takes in to account that in M33 the stars make up only ~ 57% of the baryonic mass, the 
rest being mostly in neutral and molecular gas (Corbelli 2003). 

Models of the halo occupation distribution of galaxies predict a relation for M*/M 2 oo 
as function of halo mass M200, when matching observed galaxy properties from the Sloan 
Digitial Sky Survey to the properties of dark matter halos seen in simulations (e.g., Wang et 
al. 2006; Guo et al. 2010). From SectionOwe have M M31jVir = (1.50±0.38) x 1O 12 M , which 
corresponds to Mm3i,2oo = (1.26±0.32) x 1O 12 M . Combined with knowledge of the observed 
-^M33,*/-^M3i,*, this can be used to estimate M 20 o for M33, and hence the virial mass0 This 
yields MM33, v i r = (0.170 ± 0.059) x 1O 12 M based on the Guo et al. relations, and MM33, v ir = 
(0.127±0.055) x 1O 12 M based on the Wang et al. relations. The uncertainties were estimated 
using a simple Monte-Carlo scheme that includes, in addition to the observational errors, the 
Gaussian cosmic scatter of ~ 0.2 dex in stellar mass at fixed halo mass (Guo et al. 2010)@ 
The difference in normalization between the predictions from Wang et al. and Guo et al. is 
not well understood. So we treat this as an additional model uncertainty, and allow all values 
bracketed between the two relations with equal probability. This yields as our final estimate 
M M3 3,vir = (0.148 ± 0.058) x 10 12 M @ 



7 In relating A/200 to M v i r for M33, we assume that c vlr = 10 ± 2, as we did for the MW and M33. While 
the lower mass of M33 would in principle lead one to expect a higher concentration, this is not supported 
by fits to the rotation curve (Corbelli 2003). 

8 This exceeds the observational errors in M* for both M33 and M31. As a result, it is not necessary for 
the present method to have particularly robust estimates of these observational uncertainties. 

9 We could instead have used the mass Mm33,* directly to estimate -MM33, V ir, with no reference to M31. 
This yields M M 33,vir = (0.225 ± 0.055) x 10 12 M Q based on the Guo et al. relations, and M M33)V i r = (0.123 ± 
0.034) x 10 12 Mq based on the Wang et al. relations. Both of these estimates are consistent with what we 
use here. However, there is a significant difference in absolute mass normalization between the theoretical 
relations. The relations agree better in a relative sense, which is why we prefer the method used here. The 
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Strictly speaking, the mass inferred from halo occupation distributions is the so-called 
infall mass. Thus we assume that mass loss to M31 has not yet been significant. On the 
other hand, the uncertainty on M M 33, v ir is significant. Also, our M M 33 iV ir estimate falls below 
what is implied by direct application of the Guo et al. relations. Therefore, significant mass 
loss would not be inconsistent with the range of masses we explore here. 



6.3. Mass Implications of a Bound M31-M33 Pair 

To assess the likelihood, given the data, that M31 and M33 are bound, we set up mass 
and velocity combinations in Monte-Carlo sense. The initial masses M vir for both galaxies 
were drawn as in Sections 15.21 and 16.21 The initial phase-space coordinates were drawn 
as in Sections @] and 16.11 This scheme propagates all observational distance and velocity 
uncertainties and their correlations, including those for the SuiJ^l For each set of initial 
conditions we calculated the binding energy of the M33-M31 system. The M33-M31 system 
was found to be bound in 95.3% of cases. Therefore, our observational knowledge of the 
masses, velocities, and distances of these two galaxies indicates that indeed, they most likely 
form a bound pair. 

The observation of tidal features associated with M33 independently implies that M31 
and M33 are likely a bound pair (Braun & Thilker 2004; McConnachie et al. 2009). If 
we enforce this as a prior assumption, then this affects our posterior estimates of the M31 
and M33 masses. To enforce this assumption, we merely need to remove from our Monte- 
Carlo scheme those initial conditions in which M33 and M31 are not bound. Figure Hb 
shows as solid histograms the posterior distributions after application of this additional prior. 
The main effect is to disallow some of the initial conditions in which Mm3i (blue curve) is 
on the low end of its probability distribution. The average and RMS mass increase from 
M M3 i,vir = (1.51 ±0.42) x 1O 12 M to (1.54 ±0.39) x 10 12 M o . The mass distribution of M33 
is not appreciably affected. The posterior mass distribution for M tot vir = M MWjV ir + ^M3i,vir 
is shown as the cyan histogram. Its average and RMS are 

M tot , vir = (3.17 ± 0.56) x 10 12 M (final estimate). (10) 

This is similar to the result from Section |5~2| which was M tot vil . = (3.14 ± 0.58) x 10 12 M . 
Hence, the assumption that M33 must be bound to M31 does not help much to reduce the 



latter uses only relative theoretical predictions, combined with the kinematically determined virial mass for 
M31. 

10 Uncertainties in the RA and DEC of M31 and M33 are negligible and were ignored. 
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uncertainties in the LG mass, since the fact that they are bound is already implied at high 
confidence by the observed velocities. 

The probability distributions of M31 and M33 distances and velocities are not appre- 
ciably affected by the additional prior that M31 and M33 be bound. The average positions 
and velocities in the Galactocentric rest frame remain the same to within 1 kpc and a few 
km/s, respectively, after the unbound orbits are removed. 

7. Discussion and Conclusions 

We have presented the most accurate estimate to date of the transverse motion of M31 
with respect to the Sun. This estimate was made possible by the first PM measurements for 
M31, made using HST, and presented in Paper I. We have combined these measurements 
with other insights to constrain the transverse motion of M31 with respect to the MW. We 
have used the resulting motion to improve our understanding of the mass of the Local Group, 
and its dominant galaxies M31 and the MW. 

The HST PM measurements from Paper I pertain to three fields in M31. The PM for 
each field contains contributions from three components: the M31 COM motion, the known 
viewing perspective, and the internal kinematics of M31. To correct for the contributions 
from internal kinematics, we have constructed detailed iV-body models. The models include 
both the equilibrium disk, bulge, spheroid, and dark-halo components, as well as the material 
from a tidally disrupted satellite galaxy that is responsible for the GSS. Even though the 
stars in M31 move at velocities of hundreds of km/s, the internal-kinematics corrections 
to the observed PMs averaged over all fields is quite small (< 25 kms" 1 , well below the 
random uncertainties in the measurements). This is largely due to the known properties of 
the carefully chosen field locations, and the galaxy components that they sample. 

The resulting M31 transverse motion should be largely free from systematic errors, 
based on the many internal consistency checks built into our PM program, as discussed in 
Paper I. This includes the fact that the observations for the three different fields, including 
observations with different instruments at different times, all yield statistically consistent 
estimates for the M31 COM motion. Nonetheless, an entirely independent check on the 
results is obtained by comparison to the M31 transverse motion estimates implied by the 
methods from vdMG08, which are based exclusively on the kinematics of the satellite galaxies 
of M31 and Local Group. 

Instead of using the published results from vdMG08 directly, we have redone their anal- 
ysis using expanded satellite samples, including new data that has become available in recent 
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years. The end result is similar to what was already published by vdMG08. More impor- 
tantly, the result is statistically consistent with that obtained from the HST PM program. 
Since the methods employed are totally different, and have very different scopes for possible 
systematic errors, this is very successful agreement. This gives added confidence in both re- 
sults, and also suggests that a further reduction in the uncertainties can be obtained by taking 
the weighted average of both methods. This yields (vw,vn) = (—125.2 ±30.8, —73.8 ±28.4), 
which is our final estimate for the heliocentric transverse motion of M31. The uncertainties 
in this result are similar to what has been obtained from VLBA observations of water masers 
in the M31 satellites M33 and IC10 (Brunthaler et al. 2005, 2007). 

To understand the motion of M31 with respect to the MW, it is necessary to correct 
for the reflex motion of the Sun. We adopted the most recent insights into the solar motion 
within the Milky Way. These imply an azimuthal motion for the Sun (the sum of the LSR 
motion and the solar peculiar velocity) of ~ 250 km/s, which is ~ 25 km/s higher than 
what has typically been used in previous studies. This implies a radial approach velocity 
of M31 with respect to the Milky Way of V radM 3i = -109.2 ± 4.4 km s" 1 , which is ~ 20 
km/s slower than what has typically been used in previous studies. The best estimate 
for the tangential velocity component is Vtan,M3i — 17.0 km s -1 , with la confidence region 
Kan,M3i < 34.3 km s -1 . Hence, the velocity of M31 is statistically consistent with a radial 
(head-on collision) orbit towards the MW at the la level. 

The new insights into the motion of M31 with respect to the MW allowed us to revise 
estimates of the Local Group timing mass, as presented most recently by vdMG08 and Li 
& White (2008). This yields M totitiming = 4.23 x 10 12 M for an assumed radial orbit, with 
a random error from observational uncertainties of 0.45 x 10 12 M . This result is ~ 20% 
lower than typically found in previous studies, due to the lower Kad,M3i used here. We 
calibrated the timing mass as in Li & White (2008) based on cosmological simulations. 
However, we selected from their galaxy pairs in the Millennium simulation only those with 
low Vtan,M3i, for consistency with the observations. This yields M to t, v ir = A^MW,vir±Mvi3i,vir = 
(4.93 ±1.63) x 1O 12 M , where the uncertainty now includes cosmic scatter (which dominates 
over random errors). 

We have presented a Bayesian statistical analysis to combine the timing mass estimate 
for M to t,vir with estimates for the individual masses of M31 and the MW obtained from 
other dynamical methods. For the individual masses we used relatively broad priors that 
encompass most values suggested in the literature. Even then, the cosmic scatter in the 
timing mass is too large to help much in constraining the mass of the Local Group. Its main 
impact is to increase by ~ 10% the mass estimates already known for the individual galaxies 
(and their sum). 
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In an attempt to further refine the M31 and Local Group mass estimates, we have 
studied the galaxy M33. Its known PM allowed us to study the relative motion between 
M31 and M33. A range of arguments suggests that the mass of M33 is ~ 10% of the M31 
mass. The masses and relative motions of M31 and M33 indicate that they are a bound 
pair at 95% confidence. Observational evidence for tidal deformation between M33 and M31 
suggests that the small 5% probability for unbound pairs, as allowed by the observational 
uncertainties, may not be physical. This makes low values for the M31 mass unlikely, and 
hence increases the expectation value for the Local Group mass, but only by ~ 1%. Our final 
estimate for the Local Group mass from all considerations is M to t,vir = (3.17±0.57) x 1O 12 M . 

The velocity vectors between M31, M33 and the MW are all closely aligned with the 
plane that contains these galaxies. Paper III presents a study of the future orbital evolu- 
tion and merging of these galaxies, using the velocities and masses derived here as starting 
conditions. 
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Barnes. The authors are grateful to T. J. Cox for contributing to the other papers in this 
series, and to the anonymous referee for useful comments and suggestions. 



A. Dark Halo Profiles, Masses, and Sizes 

Spherical infall models show that a virialized mass M vir has an average overdensity A vir 
compared to the average matter density of the Universe. The virial radius r v ; r therefore 
satisfies p vir = 3M vir /47rr;J ir = A vir f2 m p cri t, or in physical units (Besla et al. 2007) 

--^M^r^io^r (ai) 

For the cosmological parameters used here, h = 0.7 and fl m = 0.27, one has A vir = 360 
(Klypin et al. 2011). 

Dark halo density profiles in cosmological simulations are well described by an NFW 
density profile (Navarro et al. 1997), 

Pn{t) = p s x~ l {l + x)~ 2 , x = r/r s . (A2) 
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The enclosed mass is 

M N (r) = 4 Ws 3 /(x) = M vir /(x)//(c vir ), f(x) = ln(l + x) - (A3) 

where the concentration is defined as c vir = r vir /r s . The average enclosed mass density equals 

p N (r) = 3M N (r)/4nr 3 = 3p s (r s /r) 3 f(x). (A4) 

Another characteristic radius that is often used is the radius r 20 o so that the average 
enclosed density is 200 times the critical density of the Universe, p 200 
200p crit , where M 20 o is the enclosed mass. It follows from the respective definitions that 

Q = P 2 oo/Pvir = (200/A vir )fi-\ (A5) 

which yields q = 2.058 for the cosmological parameters used here. This exceeds unity, and 
therefore r 2 oo < ?"vir and M 20 o < M vir . For the NFW profile, r 20 o is the solution of the 
equation p jV (r 20 o) = ?Pjv( r vir), which implies 

, ( /(c 2 oo) y /3 (KfC . 

C 20 0/Cvir = — 77 x ( A6 ) 

\qf{Cvir)J 

where c 2 oo = ^2oo/ r s- This equation can be quickly solved numerically using fixed point 
iteration, starting from an initial guess for c 2 oo on the right hand side. The corresponding 
mass ratio is 

M 200 /M vir = /(c 200 )//(c vir ). (A7) 

As discussed in Springel et al. (2005), it is often convenient for numerical reasons to 
model dark halos with a Hernquist (1990) profile. This is what we will do in our exploration 
of the orbital evolution of the MW-M31-M33 system in Paper III. In this case the density 
profile is 

Ph{t) = (J^j y-\l + y)-\ y = r/a. (A8) 

Here M H is the total mass of the system, which is finite, unlike for the NFW profile. The 
enclosed mass is 

M H (r) = M H y 2 (l+y)- 2 . (A9) 
The Hernquist profile has the same density as the NFW profile for r — > if 

M H = 2np s a 2 r s . (A10) 

We can choose the scale radius a so that the enclosed mass of the NFW and Hernquist 
profiles is the same for some radius f, Mn{j-) = Mh{t). This implies 

a/r. = {[2/(x)]- 1 ^-(l/x)}- 1 , (All) 
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where x = r/r s . The corresponding total mass of the Hernquist profile satisfies 

M H /M vir = (a/r s ) 2 /[2f( CviT )} (A12) 

If we choose x = c 2 oo, then the NFW and Hernquist profiles have the same enclosed mass 
M 2 oo within r 2O o0- We denote by a 2 oo the corresponding value of a from equation (lAllj) . 



and by M# j2 oo the corresponding value of Mh from equation (1A12I) . If instead we choose 



x, = c v i r , then the NFW and Hernquist profiles have the same enclosed mass M vir within 
r vir . In this case we denote by a vir the corresponding value of a from equation flAllj) . and 
byM# )V ir the corresponding value of Mjj from equation (IA12I) . 

As an example, we consider a halo with c vir = 10. This yields c 2 oo = 7.4, M 200 /M vir = 
0.84, a 200 /r s = 2.01, M Ht200 /M vil = 1.36, a vir /r s = 2.09, and M H , V u/M vir = 1.46. 
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Table 1. M31 Transverse Velocity: Proper Motions and Internal Kinematics 







Spheroid Field 


Disk Field 


Stream Field 


HST PM Measurements 


v w (HST) 


km/s 


-167.2 ±60.2 


-194.6 ±89.8 


-65.3 ±101.5 


v N (HST) 


km/s 


-137.2 ±56.2 


-38.0 ±89.1 


-130.3 ±99.3 


M31 Internal Kinematics Model 


f (base) 




0.738 


0.922 


0.195 


/ (GSS) 




0.262 


0.078 


0.805 


vlos (base) 


km/s 


-1.1 


219.8 


4.8 


vlos (GSS) 


km/s 


71.9 


-82.3 


-185.0 


vlos (all) 


km/s 


18.0 


196.3 


-148.0 


vw (base) 


km/s 


11.0 


-25.9 


12.0 


v w (GSS) 


km/s 


22.7 


-82.6 


73.4 


vw (all) 


km/s 


14.1 


-30.3 


61.4 


vn (base) 


km/s 


-14.3 


-49.3 


-7.1 


vn (GSS) 


km/s 


2.5 


41.8 


157.5 


vn (all) 


km/s 


-9.9 


-42.2 


125.4 


M31 COM Motion (HST PMs + Internal Kinematics Model + Viewing Perspective) 


v w (COM) 


km/s 


-179.1 ±64.1 


-158.0 ±92.4 


-126.3 ±103.6 


vn (COM) 


km/s 


-122.6 ±60.0 


-0.5 ±91.3 


-247.5 ± 102.1 



Note. — Kinematical quantities for the three HST fields observed in Paper I. The top part of the 
table lists the HST PM measurements, transformed to km/s using D = 770 ± 40 kpc. The middle 
part gives the predictions of the M31 internal kinematics model described in Section [2j Predictions 
are split into two components, the "base" equilibrium galaxy model (disk, bulge, and halo), and 
the accreted "GSS" tidal stream component. The fraction / gives the amount contributed by 
each component. Predictions for all model stars (independent of component), are also listed (i.e., 
averages suitably weighted by the corresponding fractions /). Average velocities are given in the 
line-of-sight (LOS), West, and North directions. The directions are defined at the center of each 
field. The internal kinematics velocities are expressed in a frame in which M31 is at rest. The bot- 
tom part of the table gives the estimates for the center-of-mass (COM) motion of M31 that result 
from correcting the HST measurements for both internal kinematics and the viewing perspective. 
For the COM velocities, the West and North directions are always defined at the COM. 
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Table 2. Addition to vdMG08 M31 Satellite Galaxy Sample 



Name 


Type 


P 


<I> 


v\ os 






deg 


deg 


km/s 


(1) 


(2) 


(3) 


(4) 


(5) 


B517 


GC 


3.29 


77.48 


-272 ± 54 (G07) 


Mac-GCl 


GC 


3.39 


-115.38 


-219 ± 15 (G07) 


B514 


GC 


4.04 


-145.58 


-456 ± 23 (G05) 


EC4 


GC 


4.39 


135.88 


-288 ± 2 (C09) 


B516 


GC 


4.76 


28.44 


-181 ± 5 (G07) 


B518 


GC 


5.74 


-110.08 


-200 ± 48 (G07) 


And XV 


dSph 


6.84 


114.86 


-323 ± 1 (T12) 


B519 


GC 


7.35 


165.67 


-268 ± 47 (G07) 


And XI 


dSph 


7.50 


174.23 


-462 ± 4 (T12) 


And XIII 


dSph 


8.46 


166.90 


-185 ± 2 (T12) 


Mar-GCl 


GC 


8.50 


168.61 


-312 ± 17 (G07) 


And XXI 


dSph 


9.00 


-78.35 


-361 ± 6 (T12) 


And XVI 


dSph 


9.50 


158.04 


-367 ± 3 (T12) 


And XXII 


dSph 


16.06 


141.60 


-127 ± 3 (T12) 



Note. - - The sample of additional M31 satellites, which was combined with the sample from 
Table 1 of vdMG08 for the modeling of Section 13.11 Column (1) lists the name of the object and 
column (2) its type. Objects labeled "GC" are distant globular clusters. Columns (3) and (4) 
define the position on the sky: p is the angular distance from M31 and $ is the position angle with 
respect to M31 measured from North over East, calculated from the sky positions (RA,DEC) as in 
van der Marel et al. (2002). The satellites in the table are sorted by their value of p. Column (5) 
lists the heliocentric line-of-sight velocity and its error. The source of the measurement is listed 
in parentheses — Collins et al. (2009): C09; Tollerud et al. (2012): T12; Galleti et al. (2005): 
G05; Galleti et al. (2007): G07. Sky positions were obtained from the listed sources or the NASA 
Extragalactic Database. 
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Table 3. M31 Center-of-Mass Heliocentric Velocity Estimates 



Method 


VLOS 


v w 


v N 




km/s 


km/s 


km/s 


(1) 


(2) 


(3) 


(4) 


HST PMs + Internal Kinematics Model + Viewing Perspective 


(Section^) 


Spheroid Field 




-179.1 ±64.1 


-122.6 ±60.0 


Disk Field 




-158.0 ±92.4 


-0.5 ±91.3 


Stream Field 




-126.3 ± 103.6 


-247.5 ± 102.1 


Weighted Av. 




-162.8 ±47.0 


-117.2 ±45.0 


Analysis of Satellite 


LOS Kinematics (Section^) 




M31 Satellites 


-279.3 ± 16.4 


-176.1 ± 144.1 


8.4 ± 85.4 


M33 PM 


-183.1 ± 84.9 


-47.7 ± 88.2 


70.9 ± 91.5 


IC 10 PM 


-346.1 ± 84.8 


-16.2 ± 88.0 


-47.3 ± 89.3 


Outer LG Galaxies 


-361.3 ± 83.6 


-140.5 ± 58.0 


-102.6 ± 52.5 


Weighted Av. 


-281.1 ± 15.6 


-97.0 ± 40.7 


-45.1 ± 36.6 


All methods combined 


Weighted Av. 




-125.2 ± 30.8 


-73.8 ± 28.4 



Note. - Estimates of the heliocentric velocity of the M31 COM from different methods, as 
indicated in column (1). The top part of the table gives the results from the HST PM measurements, 
corrected for internal kinematics as described in Section [2j The weighted average of the results 
from the three different HST fields is listed as well. The middle part of the table gives the results 
from the updated vdMG08 analysis, based on the kinematics of M31 and LG satellite galaxies, as 
described in Section [3l The weighted average of the four independent estimates is listed as well. 
Column (2) lists the estimated M31 systemic line-of-sight velocities (the actual velocity measured 
directly from M31 itself is known to be —301 ± 1 km/s; vdMG08). Columns (3) and (4) list the 
estimated M31 transverse velocities in the West and North directions, respectively. The bottom 
line of the table lists the weighted average of the two weighted averages from the different methods. 
This is the final result used in the remainder of our study. 



3 2 10-12-? 

(degrees) 



Fig. 1. — Smoothed projected view of the N-body model used to calculate the internal M31 
PM kinematics for our HST fields from Paper I. A standard sky projection is used, with 
North up and East to the left. The GSS is visible South-East of the galaxy center, and the 
observed positions of the North-East and Western shelf are shown with dashed outlines. This 
image can be compared to star count maps of giant stars in M31, such as that reproduced 
in Figure 1 of Paper I, which show very similar features. 



-36 - 




■3oa i ■ i ■ f t : 

Projected radius (degrees) 

Fig. 2. — Smoothed view in LOS velocity vs. projected distance space of the N-body model 
used to calculate the internal M31 PM kinematics for our HST fields from Paper I. Only 
particles with Y < —0.75° are shown (located South-East of the galaxy center), where Y is 
a cartesian coordinate along the projected galaxy minor axis. The dark band in the figure is 
due to the GSS, while the base galaxy contributes most of the remaining particles. The GSS 
location matches the observed peak LOS velocity of the GSS as a function of radius (blue 
points; Ibata et al. 2004; Guhathakurta et al. 2006; Kalirai et al. 2006a; Gilbert et al. 2009), 
including that measured in the HST stream PM field of Paper I (circle). 
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Fig. 3. — Estimates of the M31 heliocentric transverse velocity in the West and North 
directions. Data points with error bars are from Table El Red: Weighted average of HST 
proper- motion measurements, corrected for internal kinematics (Section [2]). Blue: Weighted 
average of methods based on satellite kinematics (update of vdMG08 result; Section [3]). 
Black: Overall weighted average of all measurements. The starred symbol indicates the 
transverse velocity that corresponds to a radial orbit for M31 with respect to the Milky 
Way. The measurements are consistent with a radial orbit. 
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Fig. 4. — Probability distributions for the mass of the Milky Way, M31, and their sum 
M tot . The top panel shows prior probability distributions based on several lines of evidence. 
M31 (blue) and MW (green): probability distributions based on studies of these galaxies as 
discussed in the text; Sum of M31 and MW (magenta); Timing argument M to t with inclusion 
of observational errors (black), and with additional inclusion of cosmic variance from Li & 
White (2008; red). The bottom panel shows with the same color coding posterior probability 
distributions, obtained by combining constraints as described in the text. Dotted: knowledge 
of the individual galaxies combined with the the timing argument; the timing argument does 
not help much to constrain the masses, due to its large cosmic variance. Solid: Requiring 
also that M33 and M31 be a bound pair; this reduces the probability of low M31 masses. 



